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Abstract. In this paper we develop a robust multigrid preconditioned Krylov subspace method 
for the solution of heterogeneous indefinite Helmholtz problems. The preconditioning operator is 
constructed by discretizing the original Helmholtz equation on a complex stretched grid. As this 
preconditioning operator has the same wavenumber as the original problem, and the only difference 
stems from the complex stretching of the discretization grid, its spectrum closely approximates the 
spectrum of the original Helmholtz operator, resulting in a fast converging Krylov subspace method. 
We have analyzed a multigrid cycle based on polynomial smoothers and have found a condition on 
the parameters used in the complex stretching of the discretization grid such that the existence of a 
stable third order polynomial smoother is guaranteed. In practice we use three iterations of GMRES 
as a smoother, and have found this to be a viable smoother for the multigrid preconditioning process. 
We apply the method to various test problems and report on the observed convergence rate. 



1. Introduction. The efficient numerical solution of the Helmholtz equation 

Hu = ~{A + k'^)u = f in n (1.1) 

with a space dependent wave number k is important for many areas of science and 
engineering. Examples are seismic exploration |2], imaging with tomography |12j or 
photo lithography where the propagation of an optical wave through a photo resist 
determines the exposed pattern ^ . Even the Schrodinger equation that describes 
the dynamics of electrons in molecules can be rewritten as a Helmholtz equation with 
a space dependent wave number 

After discretization of 7? on a grid with grid distance h, the discrete negative 



Laplacian, — A^, in ( 1.1 1 leads to a spectrum that is spread between and a point on 
the real axis whose magnitude is of order However, the wave number will shift 

this spectrum in the negative direction. This leads to an indefinite matrix for wave 
numbers k larger than a certain threshold. This means that the spectrum is spread 
over both the positive real part and the negative real part of the complex plane, not 
necessarily excluding zero as an eigenvalue. 

This indefiniteness prevents the efiicient solution through a preconditioned Krylov 
subspace method. Typically, in these methods an unfavorable condition number 
caused by the Laplacian is improved by an efficient preconditioner. For positive 
definite problems the method of choice for preconditioning is multigrid which is linear 
in time. However, its nice convergence properties are destroyed when it is applied to 
the indefinite Helmholtz equation. 

Alternative preconditioners for the Helmholtz equation have been proposed by 
Haber and MacLachlan [TT] that transforms the problem into a complex valued 
reaction-diffusion equation. BoUhofer, Grote and Schenk developed a multilevel 
strategy based on incomplete factorizations j4|. In his PhD thesis Pinel looked at 
two-level preconditioners [13] . Gander and Nataf looked at a preconditioner which 
had a analytic incomplete factorization [TD] and Plessix and Mulder at separation of 
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variables |14j . A preconditioner based on multigrid was analyzed by Elman, Ernst 
and O'Leary [7], while Erlangga, Vuik and Oosterlee [H] used multigrid applied to a 
Helmholtz equation with a complex wave number. This paper is based on the work of 
the latter two papers. 

Two difficulties emerge when multigrid is applied to the Helmholtz problem. 
First, typical smoothers like weighted Jacobi or Gauss-Seidel are unstable for indefinite 
problems. The modes with a negative eigenvalue grow under the action of the smoother. 
These unstable modes are the smoothest modes and for slightly indefinite problems 
this instability can be compensated by coarse grid correction jSlUT]. However, for 
substantially indefinite problems the coarse grid operator itself becomes divergent and 
has been analyzed in detail by Elman, Ernst and O'Leary [7]. They found that during 
the coarse grid correction certain modes are amplified when an eigenvalue of the coarse 
grid is close to zero. This gives rise to peaks in the convergence plots as a function of 
the wave number k, which can be considered as wave number dependent resonances. 

An important innovation that avoids these resonances was introduced by Erlangga, 
Vuik and Oosterlee [Q^ that extended the work that Bayliss, Goldstein and Turkel did 
in pp. They define their preconditioner as a complex shifted Laplacian (CSL) 

M^^^u= -{A + I3^k'^)u inn, (1.2) 



where the wave number is multiplied by a complex number /3 g C. Note in particular 
that the symbol used for this complex number in ^ is /3, whereas in this paper we 
write which is more compliant with our development. For problems having a 
constant wavenumber k, the spectrum of the discretized complex shifted Laplacian 
preconditioner can be seen as a translation of the spectrum of the Helmholtz matrix 
Hfi in the complex plane. The real part of the shift defines the horizontal displacement 
while the imaginary part defines the vertical displacement. The spectrum can therefore 
be moved as far from the origin as necessary. This makes the preconditioning operator 
feasible for multigrid inversion, since the shift mitigates the coarse grid correction by 
keeping the spectrum of the coarse grid operator away from the origin. The resonances 
are now broadened and washed out. However, the shift /3 is a parameter in the 
method. From an iterative viewpoint the preconditioner and the original problem are 
inversely related. Measures to enhance the multigrid inversion (by increasing the shift 
magnitude) deteriorates the Krylov convergence, and vice versa. The optimal shift is 
therefore a compromise between the multigrid and the Krylov parts of the solver. 

The eigenvalues of the resulting discrete preconditioned system, (M^^^)~-^ Hh, lie 
inside a small circle in the complex plane, away from zero, when Sommerfeld radiation 
boundary conditions are applied. The same complex shift idea can be used for other 
absorbing boundary conditions such as exterior complex stretching (ECS), where the 



Helmholtz operator in ( 1 . 1 1 is defined on a partially complex domain |15j . Similar to 
Sommerfeld conditions the spectrum of the preconditioner is then shifted away from 
the origin, while the preconditioned system is ideally clustered for the outer Krylov 
subspace method. In [TS] the alternative complex stretched grid (CSG) preconditioner 

M'^^'^u = -{A + k'^)u inn (1.3) 



is defined, where the Helmholtz operator is retained but the domain is altered so that 
A is discretized with a complex valued grid distance. A proper choice of the domain 
n, typically a complex rotation over an angle introduces a favorable rotation in 
the spectrum of the Helmholtz operator, away from the origin. 
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The introduction of the complex shifts, either in (1.2) or (1.3 1, avoids the ap- 
pearance of resonances but it does not prevent the smoother from becoming unstable. 
In addition, both preconditioners leave the user with choosing a proper shift /3 or 
domain rotation 6p, respectively. This papers focuses on finding a suitable smoother 
for complex shifted or stretched Laplacians with ECS absorbing boundary layers. We 
focus on polynomial smoothers that can be viewed as products of multiple w-Jacobi 
steps with a different choice of lo in each step. Requiring that this smoother is stable 
leads us to conditions on the parameters used in the preconditioner. 

In the numerical results section we report also the convergence rates when the 
polynomial smoother is replaced by GMRES(m). And we find that using GMRES(3) 
on all levels gives a satisfactory convergence. This builds on the work of Elman et al. 
|71 which analyzed and optimized a smoother schedule that combines Jacobi smoothing 
with GMRES. However, this gives rise to complicated multilevel methods, where many 
GMRES iterations are required on some levels while on others a Jacobi sweep suffices. 

In this paper the wave number is denoted by k. M^^'^ denotes the discretized 
complex stretched grid preconditioning matrix, while M^^'^ denotes the continuous 
preconditioning operator. Furthermore we use the following notations, see also Fig- 
ure h for the original interior real grid distance of the ECS grid, 9^ for the rotation 
angle of the exterior absorbing ECS layer, the exterior complex grid distance of the 
ECS layer, = he^^i and Op for the rotation angle of the interior complex stretched 
domain for M'~^^'^ . In the same way we use hp the interior complex grid distance 
for Mj^^^, hp — he^^f. We also use A^'^-* for the eigenvalue corresponding to the 

eigenmode with the highest energy, or most oscillatory, on the interior region, A^"^ 
for the eigenvalue corresponding to the eigenmode with the highest energy, or most 
oscillatory, on the ECS layers and Aq ~ — A:^ the lowest energy eigenvalue, also referred 
to as the smoothest eigenvalue. 
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Fig. 1.1: A one-dimensional illustration of the computational grids, the extension to 
multiple dimension is straightforward with Kronecker products. The original ECS grid 
(o) for the Helmholtz operator H has a real grid distance h in the interior region, for 
< a; < 1, and a complex grid distance on the ECS layer, for a; > 1. The CSG 
domain (+) for the complex stretched grid preconditioner M^^'^ has a complex grid 
distance hp in the interior region, for < a; < 1, and the same complex grid distance 
h-y on the ECS layer, for a; > 1. 



2. Background on the spectrum. In this section we will briefly summarize 
the results from [T^ for a D-dimensional model Helmholtz problem. In the next 
section a lower bound for the CSG-parameter 9p is found in order to have a stable 
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multigrid method for the preconditioner resulting in a good overall convergence of the 
preconditioned Krylov subspace method. 

Suppose we have a one dimensional complex grid 



_ Jh (0<j<n), 

I 1 + — n)n^i [n+1 < ] <n + m), 

that consists of n intervals of complex grid distance /i^ followed by m intervals of 
complex grid distance as illustrated in Figure |1.1[ 



This particular kind of grid is used to create a CSG preconditioner (1.3) by 
discretizing the Helmholtz equation on it. It originates from an ECS grid for the 
original Helmholtz matrix Hf^ by replacing its interior real grid distance hhy hp = he^^f 
for complex rotation or hp = h(l + itaii(6p)) for complex shifting [15 . The exterior 
grid distance hj — he^^'' or h^ = h{l + ii'An{0^)) is the same for Hh and M^^'^ and 
forms the absorbing ECS layer. 



We discretize the Laplacian operator A in Equation (1.1 1 on (2.1) with the 



Shortley-Weller finite difference scheme for non-uniform grids in one dimension 
(Pu 2 f I / 1 1 \ 1 



dz"^ hj^i + hj \hj^i V^j-i J ^ 



[zj) ~ , , , 1—^,-1 - ^ + ^ «J + (2-2) 



in grid point j, where hj^i and hj are the left and right mesh widths respectively, and 
may belong either to the hp category or to the h^ category. The discretization leads 
to a linear system 

M^^^Uh = -{Lh + k'lh)uh = 9h, (2.3) 

where Lh represents the discretized Laplacian, 1^ is the identity matrix and gh is a 
vector containing source function values sampled at the grid points. The discretization 
matrix for the higher dimensional problem is created by Kronecker products of the 
one-dimensional matrices. 

Lemma 2.1. The spectrum of the D-dimensional Helmholtz operator in Equa- 



tion (1.1) discretized with the Shortley-Weller formula (2.2) on the CSG-grid defined 
in (2.1), is hounded in the complex plane by a triangle zqZ\Z2 described by the complex 
points zq = — fc^, zi = —k^ -\- 4,D/h'p and Z2 = —k^ -\- 4Z)//i^. 

Proof. This easily follows from the results in [T5] for the discretized Helmholtz 
equation on an ECS-grid. □ 

The eigenvalues of the Helmholtz discretization matrix on the CSG-grid are not 



randomly distributed inside the triangle zoZiZ2, defined in Lemma 2.1 their exact 
location is related to a physical interpretation. The eigenvalues that are close to 
the vertex zq = —k^ are aligned and correspond to the smoothest or low frequency 
eigenvectors on the entire CSG-grid. We will refer to them as the smoothest eigenvalues. 
There are oscillatory or high frequency eigenvectors localized on the interior part of the 
CSG-grid with eigenvalues in the neighborhood of the second vertex zi = —k'^-\-AD/h'p. 
The region around the third vertex Z2 = —k"^ -\- iD/h^ contains the eigenvalues with 
high frequency eigenvectors localized on the absorbing ECS layer. We will refer to 
them as the high frequency eigenvalues on the interior part of the domain and the 
ECS layer respectively. 

Remark 1. We will denote the three eigenvalues that are closest to the vertices by 
'^a ~ Zq, Ajif"* ~ zi and X'^^ ~ Z2. They move closer to the three vertices of the triangle 
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with finer discretization, so for simplicity we will no longer make the distinction 
between the vertices zq, zi and Z2 and the actual extreme eigenvalues Xq, A^'^'' and A^y^-* 
in the remainder of the paper. 

3. A Polynomial smoother. Building a multigrid method involves designing 
a hierarchy of grids of different resolutions together with the proper restriction and 
prolongation operators to translate information from one grid to another. Another 
important ingredient of the multigrid method is the smoother. This is the relaxation 
method that is used on each level to eliminate the high frequency modes of the error. 
Multigrid methods can be very efficient when smoothing and coarse grid correction 
complement each other and perform poorly (if at all) when this complementarity 
cannot be put together. Due to indefiniteness the discretized Helmholtz problem poses 
this difficulty for multigrid built on standard relaxation methods such as w-Jacobi or 
Gauss-Seidel. Even if the discretization matrix has all nonzero eigenvalues on the finest 
multigrid level, one or more coarser representations can have eigenvalues undesirably 
close to zero that can destroy the smoothing of the relaxation method on that level and 
therefore the entire multigrid performance. One workaround is to invest computational 
effort in more robust smoothers. In [7] classical smoothers were replaced by the 
GMRES method on problematic levels when multigrid is applied to the Helmholtz 



problem (1.1 1. Instead of running a few iterates of a relaxation method, GMRES is 
applied to the error equation. When this is applied to the original Helmholtz problem 
it can take a large number of expensive iterations on every troublesome level. In the 
end this GMRES accelerated multigrid method is used as a preconditioner for an outer 
GMRES method on the Helmholtz problem. 

In this paper we explore the idea of the GMRES accelerated multigrid method but 
applied to an altered Helmholtz operator such as the preconditioner M'-^^'^ , instead of 
to the original H. The advantage is that it needs at most three GMRES iterations 
to cope with potentially indefinite multigrid levels. In this section we defend this 
statement with theoretical results for the Mf^^'^ preconditioning matrix by introducing 
an idealized polynomial smoother. 

3.1. Construction. We will construct a third order polynomial smoother f{t) G 
PsiC) with complex coefficients. It is intended to work for the Helmholtz problem 



discretized on a complex stretched grid as in (2.1 ), which defines the preconditioning 
matrix Ad^^'^ . This polynomial should have the following properties to fit the multigrid 
philosophy. 



First, it should be stable: VA e cr(Mp<5) : |/(A)| < i 



• Next, it should smooth the most oscillatory modes efficiently. The mode with 
the highest frequency in the interior region has an eigenvalue A^^'^^ = — fc^, 
for a _D-dimensional problem. The highest frequency mode on the complex 
absorbing layer has an eigenvalue A^'^'' — — k^. So we should demand 

that f{t) maps these eigenvalues as /(A^'^'') = and f{X^^) — 0. 

• In contrast the smoother should leave the smoothest mode, the eigenvector 
that has an eigenvalue A — —k'^, untouched. So it is required that this mode 
is mapped by the polynomial to the unit circle, i.e. /(— fc^) = e"''. 

• Finally, it is natural that the solution of the linear system is a fixed point of 
the smoother, so /(O) = 1. 

It will be shown that such a polynomial smoother can be constructed for each level in the 
multigrid hierarchy if the complex stretched grid for Mp^'^ fits certain requirements. 



The third order polynomial can be written as 



fit) = (1 - - C^2t)(l - ^2t), (3.1) 

where uji, uj2, W3 G C. This polynomial can be interpreted as a sequence of three 



w-Jacobi steps with different weights. The coefficients uj2 and uj^ in (3.1 1 are set by 
the requirements thai 
the eigenvalues associf 
zeros of f{t) we have 



the requirements that the most oscillatory modes should be mapped to zero. Since 

X^*^^ and A^,''' 



the eigenvalues associated to the most oscillatory eigenvectors Ai*^-* and A^,'"'' should be 



/(t) = (l-^it)(l-t/A('=))(i_i/AW) 



We now choose Wi such that /(— fc^) = e*"^ with (p e [0,2tt). This ensures that the 
smooth modes are hardly affected by the polynomial smoother. This leads to 



So the choice of the point on the unit circle to which the smoothest eigenvalue A = — fc^ 
is mapped by f{t) defines the coefficient uj^p. However, as we will discuss next, not 
all choices of (p lead to a stable smoother. We now treat (p as parameter of the smoother. 

Example 1. To illustrate how the spectrum is mapped by the third order polynomial 
for various choices of ip we discuss the example that is shown in Figure [gjj for a wave 
number k — 100. The operator is discretized, as discussed in section^ with 192 finite 
difference points in the interior and 64 points on the exterior in each dimension. This 
gives a grid distance h — 0.7812 • 10~^ that is rotated in the interior by 8 degrees to 
hp — e**Tso/i and in the exterior by 30 degrees to h^e^'^^wh. 

The largest eigenvalues of M^^^ are x''^ = (11.5 — z3.61) • 10^ and xi^^ = 
(5.55 — ill. 35) • 10^. The spectrum lies in the third and the fourth quadrant of the 
complex plane and is bounded by the three line segments that connect 



two eigenvalues. This is shown in Figure 3.1 a) 



In Figure\3l\c), we show the results for p = 0.015, 0.016, 0.017, 0.018, 0.019, 
0.020, 0.021 radians. For each of these angles we have different function /, only for 
larger (/? the map is stable. 

(k) (k) 

Remark 2. We have made the choice to map X^ and X-y to zero. This does not 
necessarily lead to the most efficient smoother. It might be better to choose eigenvalues 
that are slight inside the triangle to map to zero. In a similar way for Lo-Jacobi a 
choice of uj — 2/^ leads to a better smoother than a choice of uj ~ 1/2 for the Poisson 
problem with Dirichlet conditions J^. 



3.1.1. Explicit choice for Lp. This section shows that on all levels of the multi- 



grid hierarchy a choice for ip in (3.2 1 can be made such that the smooth modes are 
mapped inside the unit circle. We represent the polynomial in (3.11 in the general 
form, 



f{t) =a + b{k^ + t)+ c{k^ + tf + d(fc2 + tf 



(3.3) 



but the coefficients are, at this moment, unspecified. The polynomial can be viewed 
as a Taylor expansion around — fc^ and so we have a simple relation between the 
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coefficients and the derivatives 



a=/(-fc2), 

c=/"(-fc2)/2, 
rf=/"'(-fc')/6. 



We will find a, b, c and d e C so that the following requirements are met: 



/(O) = 1, (fixed point property) 
f{—k'^) = e'"^ , (map —k^ to the unit circle) 
arg(/'(— fc^)//i|) = <y9 — 7r/2, (map of top edge is tangent to the circle) 

(high energy eigenvalues mapped to zero) 



/(A^''-') = 0, (high energy eigenvalues mapped to zero). 



(3.4) 
(3.5) 
(3.6) 
(3.7) 
(3.8) 



These conditions rule that the exact solution should be a fixed point (3.4 1, and that 



the smoothest modes should be mapped to the unit circle (3.51. In (3.6 1 we make 



the choice to select (p such that the top edge of the bounding triangle, that can be 
parameterized by — fc^ + s//i^ with s € [0,8], is mapped tangent to the unit circle in 
e'"^ by f{t). Equations (3.7) and (3.8) rule that the most oscillatory modes should be 
mapped to zero. This elaborates upon the idea of an ideal smoother: stability for all 
eigenmodes and the most error reduction in the range of the most oscillatory modes. 



We immediately see that a 



An alternative way to express that the first 



order term of / in —k should touch the circle is to impose that a and bs/hp £ C 
should be orthogonal when we translate to K^. This leads to 

^{b/h}) cos((^) + '^{b/h}) sin(^) = 0. (3.9) 

Since both b and f are unknowns this is a non-linear equation. The solution will give 
bounds on (p. 



Lemma 3.1. Let k be the wave number of the D- dimensional Helmholtz equation 
and hp and h^ be the mesh widths on the interior and the exterior region, respectively. 
Let (6i,ci,di)"^ and (&2jC2,d2)"^ be the solution the linear system 

e fc4 ^6 N 

Ad/hl {MY /hi {Mf/hl \x = y (3.10) 
4d//i2 (4d)V/i^ (4d)3//i^ / 

with right hand sides y = (1, 0, 0)^ and y — (1, 1, 1)"^, respectively. 
If 

mb2/hi)\<\b,/hi\ (3.11) 

then there are two solutions for ip, 

V? = arg(6i/ft2)±arccos(^l^^). (3.12) 
7 



Proof. The system of equations represented by Equations (3.4)-(3.8l has five 
equations for the unknown coefficients a,b,c,d and (p. Equation (3.5 1 immediately 
impHes a = e"^, where the parameter (p e [0,27r) is still to be determined. The 
following strategy solves for the remaining coefficients. We take the first and the last 
two equations and solve the resulting linear system 



/ 



fc4 



fc6 



fc2 + A 

fc2 + A. 



(fc) 

p 

(fe) 











1 

















(3.13) 



Note that fc^ + A^''^ 



4L>//i2 and fc^ + A^''^ 



AD/hl 



For every angle ip this system has 

a unique solution if Xf^ ^ Xil'^ ^ 0. It is a Vandermonde matrix that appears typically 
in interpolation problems. The solution can be written as a linear combination of the 
solutions (6i,ci,di)-^ and (62,22,^2)^ of the same linear system with (1,0,0)"^ and 
(1, 1, 1)"^ as right hand side. The solution for b, the first derivative of /(<) in — fc^, is 
then hi - e"^b2. 

With this solution we can find an angle </? that satisfies the remaining condition 



(3.61, 



^{b/hl) cos(^) + sin((^) = 0, 

^ - e''n2)lhl) cos(</7) + 5((5i - e'^b2)/hl) sin(v5) = 0, 

^ ?ii{bi/h}) cos((p) + sm{ip) - ^{b2/hl) ^ 0. 

If |5ft(62//i^)| < then this last equation for (p has solutions 

ip = arg(6i//io) ± arccos( ., , ). 

I^il 





(3.14) 



(3.15) 



Lemma 3.1 gives two possible values for tp in (3.2) such that the top line segment 
of the bounding triangle of the spectrum of the preconditioner is mapped tangent to 
the unit circle. However, one choice of ip gives a stable spectrum while the other does 
not. 



Corollary 3.2. //A^''' ^ Xij'^ ^ and |3fi(62//i^)| < then only one of 

the two choices of (p guarantees that f(t) is stable for the smooth modes. 



Proof. Since A^f-* 7^ A^,'^'' ^ fc^, the system (3.4l-(3.8) has a nonzero solution. The 

s/h'i with s e [0,8], 



7- '^7 

respective maps of the two line segments, — fc^ + s//i| and —k 
of the bounding triangl e tha t meet in —k'^ also meet in f{—k^) 
for ip given by Lemma 



3.1 



The two choices 

map the line segment that connects —k^ and A^'^'', which 
is closest to the real axis, such that it is tangent to the unit circle in e"^. Since the 
angle between this line and the line segment that connects —k and A^ is preserved 
by the holomorphic polynomial map /, this second line is either mapped outside the 
unit circle or inside the unit circle. The low energy eigenvalues near —k'^ are bounded 
by these two lines and so these eigenvalues are mapped inside the unit circle if this 



second line is also mapped inside the circle. □ 



Remark 3. Note that if the system (3.4l-(3.8l has no solution, i.e. |5R(62//i^)| > 



then the inner product of the first order term f^—k"^) and the second order 



term f'{-'k'^)s/h'p, (3.15) in the proof of Lemma 



3.1 



is either positive or negative. 
This corresponds to mapping the top line of the bounding triangle outside or inside 
the unit circle respectively, instead of having it as a tangent. However for all physical 



choices of the ECS-layer, and thus and X''^' , we have 3fJ(&2/^|) > and so the 
inner product will be negative for all possible choices of tp, which will always result in 
an inward mapping of the top line of the triangle. In this case we pick ip = and Wi 
is then 



Wi = 




(1 



kyxf^){i 



Example 2. It is useful to check the condition in (3.111 of Lemma 3.1 in a real 



multigrid context. In this example we look at a hierarchy of a two-dimensional problem 
where the mesh widths in each dimension are hp = 2~'e*^'^ in the inner region and 
hj = 2~'e'^^ on the absorbing layer. On each multigrid level I — 1,2, .. .7 we need to 
construct a stable smoother. 

For the fine levels, hp and h^ are small and X^p^ — — fc^ + 8/h'p and x'^'^ = 

—k^ + 8//i^ are dominated by the second term. We then have that A^*^-* « 22'+'^e^'^'^ 
and Xi^"^ ~ 22'+3g-«ST ^yj^; both lie in the fourth quadrant of the complex plane. Since 
the smooth modes lie near — fc^ in the third quadrant, the spectrum of the preconditioner 
matrix Alj^^^ on the finest levels extends from the third to the fourth quadrant. 

On the coarsest levels, however, the first terms in x'j^^ = — fc^ + 8/h'p and X^"^ = 
— fc2 + 8//i^ is small and the shift over —k^ dominates. So for the lowest levels it is 

possible that both A^'^'' and x'^^ lie in the third quadrant of the complex plane. The 
spectrum will then be negative definite as it fits completely inside the third quadrant. 

Since hp and h^ change over the levels the solutions bi and 62 of the linear system 
in Lemma 3.1 will change as well on each level. Consequently, tpi along with uj^p will 



also change from level to level. Note that these are the coefficients of the smoother in 
Equation \3.S\ 

In Figure 3.2 we check the condition in (3.11) for each of the levels. The dashed 
line marks where 3?(&2/^^)| — \bi/h'jp\ as a function of k andOp, the rotation angle of 
the mesh width hp. The dashed lines to the left come from the coarsest levels, while 
the lines to the right come from the finer levels. For a particular level I and for the 
wave numbers k on the right of the corresponding dashed line we can choose any ip 
and get a stable smoother. For wave numbers on the left of the dashed line, we have to 
solve for ip as discussed earlier in order to achieve a stable smoother. 

The other way round, for a given wave number, the ip of the coarsest levels can be 
chosen arbitrarily since the condition of Remark^is fulfilled. For those levels we pick 
Lp — Q and u)^p accordingly. This example is continued in Example [5| 

3.1.2. Stability for all eigenvalues. The condition that the image of the line 
between — fc^ and X'p^"' is tangent to the unit circle is not sufficient to guarantee stability. 
In Figure [373| we illustrate that, although, the image is tangent to the unit circle, it 



can lead to an unstably mapped spectrum. In this section we want to extend the 
region in the spectrum of smoothest eigenvalues that are mapped inside the circle. 
Again, we look at the line 

-fc^ + s/h^, with s e [0,8], 

where /i| is the mesh width of the inner part of the domain. It is mapped by the 
polynomial to 



We want this to stay inside the unit circle so we have that 

|/(-fc2 + s//i^)|2 < 1 for s<l 

Lemma 3.3. Let k be the wave number of the Helmholtz problem and hp and 
the mesh widths of the CSG-preconditioner such that there is a stable solution of the 
tangent equation for ip. Let b and c be the coefficients of the corresponding polynomial 



of the form (3.3 1. // 



'h^ 
lip 



')<-; 



hi 



then an extended part of the spectrum is mapped inside the unit circle. 
Proof. The condition |/(— fc^ + < 1 becomes 



(3.16) 





s^ + 0{s^) < 1. 

(3.17) 

for s ^ 1. Since a — e^f we have that aa* — 1. Because is such that p-a*+^p-^ a = 

and the line segment that connects —k and \\ is mapped inside the circle, it is 
clear that 




25R(-^a*) + 



hi 



< 



is required to ensure a stable map. □ 



Remark 4. // the condition in (3.161 is fulfilled, the set of eigenvalues that is 
mapped inside the unit circle is extended. However, to guarantee that the complete 
spectrum is mapped inside the unit circle, we would need to check the higher order term 
in ^J7\. Since f{-k^ 



s/hl) is zero for s — 8 and f{t) is a low order polynomial, 
it is likely that the condition in (3.161 is sufficient to make the map of the complex 



spectrum stable. Numerical experiments confirm this conjecture. 

Example 3. This example is a continuation of Example^ and it uses the same 
multigrid hierarchy to illustrate how we can use the extra condition in \3.1(^ We extend 
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the plot of Figure \372\ and add the condition of Lemma \3.3\ The solid line in the figure 

shows where 25R(^e~"^) + A- — 0, as a function of k and the rotation angle Or of 

"■p "-fi 

hp for each level in the multigrid hierarchy. In red we have added the condition of 
Lemma \37^ We see that for a given k the smoother is only fully stable for a rotation 
angle 9p above the red line. In principle we can adjust hp of the CSG until this 
condition is met. 

For example at k — A we see that the three coarsest levels have their bounding 
line to the left. That means that we can pick any angle Lp for these levels. For the 
fourth level the bounding region is on the right. This means that we have to solve 



Equation 3.12 however, if the rotation angle 9p is smaller that 5 degrees, we are not 
.sure that the spectrum is completely mapped inside the unit circle. We therefore pick 
the rotation angle of hp to be larger than 5 degrees. For the finer levels we are still 
inside the region where we need to solve for tp, but the angle of hp can be taken very 
small here. 

This example also illustrates how multigrid on the original Helmholtz matrix will 
behave badly. Indeed, we see that the extra condition on 9p for the preconditioner is 
necessary to ensure stability. The drawback is that we need to choose the same hp for 
the complete multigrid hierarchy. So in this example 9p, and thus hp, is determined 
by the fourth coarsest level. 

Remark 5. Figure \3^ illustrates that for particular choices of k of the Helmholtz 
problem the rotation angle 9p of the complex stretched grid preconditioner can be very 
small. For these cases M^^'^ and Hh will be very close to each other and the Krylov 
iteration will be efficient. However, for realistic problems when the wave number is 
space dependent it is very hard to identify these wave numbers in a robust way. For 
robustness it is therefore a good strategy not to choose the smallest possible 9p, since 
small variations in the wave number will turn the smoother unstable. 

3.2. Strategy to determine the parameters of the smoothers for all 
levels. From the examples it is clear that a stable third order polynomial smoother is 
easy to construct for all the levels except one. Indeed, let us assume that we have to 
solve a problem with fc = 75 and that we use the multigrid hierarchy of Example [2] 
and [3] For the four coarsest levels this k is on the right of the dashed line, which 
means that for these levels = is a stable choice. For these levels the spectrum of 
the CSG-preconditioner lies almost entirely in the third quadrant. In a similar way 
the smoother of the finest levels is easily constructed since for those levels the solid 
line leads to negative angles. However, for level Z = 5 we can only guarantee a stable 
smoother if 9p > 5.5. So once we have determined hp we need to calculate ip and the 



corresponding uj^ in (3.2 1 for each level. For the coarsest level we can pick any ip but 
we prefer tp — 0. 

While this strategy leads to a stable smoother for the homogeneous Helmholtz 
problem, a more robust smoother is preferred for heterogeneous problems, where 
the spectrum of the discretization matrix could slightly deviate from the describing 
triangle. Since our manually tuned third order polynomial is a particular polynomial 
with p{0) = 1 it gives an upper bound for the convergence of GMRES(3) as a stand 
alone relaxation method. As a consequence a multigrid method for a preconditioner 
such as M^^'^ only needs three iterations of GMRES when it is used as a smoother. 
This is in contrast to a much higher and undetermined number of iterations for the 
original Helmholtz matrix Hh as is shown in [7j. 

To justify the smoothing property of GMRES (3) consider again the bounding 
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triangle in Lemma [2. 1[ The lowest energy eigenvalues in the top of the triangle, around 
= 0(1), are relatively close to the origin compared to the high energy eigenvalues 
Xf) = + Ad/hj = 0{l/hl) and A^''^ = -k^ + Ad/Ki^ = 0(l//i2), especially when 
a high accuracy condition h = 0(k~^^^) is used. So whenever the minimal third order 
polynomial p{t) is constructed for the spectrum with the constraint p(0) = 1, then 
the region of the smoothest modes will be pulled towards 1 by this polynomial map 
whereas the high energy region is more freely minimized. 

4. GMRES as smoother. 

4.1. Two-grid analysis. We focus our initial numerical experiments on solving 
the preconditioning problem M^^'^Uh — bh efficiently with multigrid. First, we collect 
the convergence rates of the two grid operator 

TG = S^^il - l'^{M§^^y^I^Ml^'^^)S''\ (4.1) 

where S denotes the polynomial smoother, and the restriction and interpolation 
operators, respectively. The coarse grid operators M^^'^ is a rediscretization of the 
problem with grid distance H, as discussed in the previous sections. 

In the numerical experiments we solve M^^'^Uh = hh with zero right hand side 
and a random initial guess uf'^. The reported convergence rates are the asymptotic 
rates measured by the ratio ||e'^'^-'||/|je'-'^~^-*|j. We compare the convergence rate of the 
polynomial smoother, with the parameters chosen as in the previous section, with 
GMRES(l), GMRES(2) and GMRES(3). The measured convergence rates for various 
levels in the multigrid hierarchy are shown in the Figure |4.1[ 

Note that the explicit constructed third order polynomial forms an worst case 
bound for the convergence of GMRES (3). 

4.2. Performance of a V-cycle. After the two grid analysis multigrid with 
GMRES as smoother is tested on the preconditioner M^^'^ for the homogeneous 



Helmholtz problem ( 1.1 1 in a square domain U, = (0, 1) with a point source 



<5(.,2/) = l'' ^^^=1/2 = ^' (4.2) 
I 0, elsewhere. 

All boundaries are extended with an ECS-layer with an angle 9^ = ^ to absorb 
outgoing waves. Although a more stringent relation between k and the mesh width h 
might be desirable for physical accuracy, we only wish to observe iterative behavior 
for which 10 points per wavelength suffices, which translates into the requirement 
kh < 0.625. The preconditioning matrix M'j^^'^ is constructed by discretizing the 
problem on a complex stretched grid with a small inner angle 6 p — Q .18 k, > 
where 6*0 is the critical angle for stability of the polynomial smoother. The angle for 
the EGS-layers is kept at 6*^ = ^ as in the original Helmholtz problem. 

The preconditioning matrix M^^'-^ is inverted by applying subsequent \{v\,V2)- 
cycles where the standard smoother, e.g. w-Jacobi, is replaced by the GMRES method, 
where vi and V2 indicate the number of pre and post smoothing steps, respectively. 
One smoothing step consists of a GMRES solve on the error which is stopped after 
s iterations, simply denoted as GMRES(s). It constructs a minimal polynomial of 
order s and so will cost s matrix vector products. The number of cycles needed to 



converge to a relative residual of order 10 ^ are presented in Table 4.1 for different 
configurations and wave numbers. In the experiments with the V(l,0)-cycle we see that 
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(i^l,l'2) 


GMRES(s) 


fc 
on 




OU 


on 


CPU (fc 80) 


(1,0) 


s = 2 


18 


19 


20 


20 


7.23 




s = 3 


13 


14 


14 


14 


7.88 




s = 4 


11 


12 


12 


13 


8.58 


(1,1) 


s = 2 


11 


11 


10 


12 


7.85 




s = 3 


8 


9 


9 


9 


8.53 




s = 4 


7 


7 


7 


7 


10.33 


(2,1) 


s = 2 


9 


9 


8 


10 


17.46 




s = 3 


7 


7 


7 


8 


44.35 




s = 4 


6 


6 


6 


6 


94.13 



Table 4.1: Number of multigrid V-cycles with different number of GMRES smoothing 
steps to solve the preconditioner for a constant k. 



s — 4 does not improve the method much more than s = 3. Earher we showed that a 
third order polynomial is sufficient to guarantee a stable smoother. The experiments 
with s = 2 also converge, but only after a substantially larger number of V-cycles. 
Increasing the number of smoothing steps (vi,i>2) to more than (1,1) does not pay 
off in the eventual number of V-cycles needed to reach the tolerance. The table also 
shows the timings for k = 80. The V(1,0) with 2, 3 and 4 and V(l,l) with s = 2 and 
s = 3 perform similarly. 

Based on these observations we prefer V(1,0) and V(l,l)-cycles with GMRES (3) 
smoothing for the numerical experiments presented in the next section. Figure |4.2| 
shows that the averaged convergence rate of these V-cycles do not grow larger than 
40% (-I-) and 25% (*) for V(1,0) and V(l,l), respectively. Note that there is, however, 
a slightly variation as a function of the wave number k. The convergence rate of the 
first V(1,0) or V(l,l)-cycle alone is below 35% and 8% respectively. In the numerical 
experiments only one V-cycle will be used to approximately invert the preconditioning 
matrix M^^'^. The convergence rate of the V(2,l)-cycles is relatively close to that of 
the V(l,l)-cycles as could be expected after considering Table 



4.1 



5. Numerical experiments. In this section we test the M'~^^'~^ preconditioner 
on three two-dimensional benchmark problems. The experiments cover both homoge- 
neous media with constant wave numbers k and heterogeneous media where the wave 
number is space-dependent. The discretized problem is solved with a Krylov subspace 
method that we will call the outer Krylov method. The complex stretched grid matrix 
Mf^^'^ is constructed from the same original Helmholtz equation but on a different 
grid with the inner part slightly rotated in the complex plane over an angle Ojj. It is 
used as a preconditioner and approximately inverted with one V-cycle with GMRES (3) 
smoothing. As a consequence of this non-standard smoother the actual preconditioner 
is not the same in every outer Krylov step. Therefore the flexible GMRES method [T3] 
and the Bi-CGSTAB method are used as outer Krylov subspace methods. Whenever 
Bi-CGSTAB converged after the first matrix-vector product of the last iteration, this 
iteration is counted as a half. The experiments are all run in Matlab on two quad core 
Intel Xeon CPUs (E5462 @ 2.80GHz). The mentioned CPU times are scaled by the 
total number of grid points, including the ECS boundary layers, so it expresses the 
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computational cost per grid point. 

Case 1: Constant k.. In the first experiment we continue with the homogeneous 
Helmholtz problem with the point source in (4.2 1 from section 4.2 The complex 
shifted Laplacian preconditioner was tested on this model problem for Sommerfeld 
boundary conditions in 8 . The idea was also compared with the complex stretched 
grid preconditioner when applied with ECS boundary conditions in [15' . For the latter 
boundary conditions both strategies result in a good preconditioning matrix for a 
Krylov subspace method but still need to be (approximately) inverted by a cheap 
method. A multigrid method is preferred for this purpose as it is easily extended to 
higher dimensions where other standard methods, such as ILU, suffer from severe 
memory problems. Automatic configurations of multigrid (in this context) are often 
less efficient than manually tuned ones. In this paper we improve the multigrid 
performance on the complex stretched grid preconditioner in order to achieve a better 
overall convergence and to develop a robust and efficient solver for indefinite Helmholtz 
problems. 





k 


20 


40 


60 


80 


100 




interior grid 


32^ 


64^ 


64^ 


128^ 


256^ 


MG(V(1,0)) 


FGMRES 


26(0.08) 


46(0.7) 


54(1.65) 


77(2.78) 


84(4.03) 




Bi-CGSTAB 


19.5(0.10) 


35.5(0.61) 


50.5(1.64) 


82(2.69) 


99(3.39) 


MG(V(1,1)) 


FGMRES 


17(0.08) 


26(0.48) 


34(1.27) 


44(1.79) 


47(2.27) 




Bi-CGSTAB 


9.5(0.09) 


17.5(0.52) 


22(1.29) 


35(2.04) 


32.5(1.97) 


MG(V(2,1)) 


FGMRES 


14(0.15) 


21(2.34) 


28(6.17) 


37(8.34) 


39(5.89) 




Bi-CGSTAB 


8.5(0.17) 


15(1.86) 


19(5.04) 


27(2.33) 


26(6.03) 


ILU(O) 


GMRES 


66(0.18) 


147(3.7) 


266(20) 


294(24.1) 


490(96.1) 


ILU(O.l) 


GMRES 


77(0.20) 


160(4.2) 


274(22.3) 


317(30.8) 


1722(1309) 



Table 5.1: Number of FGMRES, Bi-CGSTAB and GMRES iterations (and CPU time 
per grid point xlO^) for different preconditioning strategies for Case 1. 



Table |5.1| displays the number of iterations and the CPU time per grid point 
xlO^ for FGMRES and Bi-CGSTAB for different values of the wave number k. The 
number of outer Krylov iterations with MG(V(1,1)) is significantly smaller than 
with MG(V(1,0)) for all values of k and results in a faster CPU time. Although a 
MG(V(2,1)) inversion of the preconditioner requires even less outer Krylov iterations, 
the total CPU time is again higher than for the MG(V(1,1)) case. This confirms what 
was observed in Table |4.1| and Figure |4.2[ that extra smoothing steps do not pay off 
because it does not quite improve the approximate preconditioner inversion. Figure [O] 
visualizes the typical relation between the wave number k and the convergence speed, 
which is obviously influenced by the increasing number of grid points for growing k. 

As a reference the same problem is solved with ILU(O) and ILU(O.l) inversion of 
the preconditioner and GMRES as the outer Krylov method. Both iteration numbers 
and CPU times are not competitive to multigrid inversion of the preconditioner. 

Case 2: Wedge problem. In the second experiment the proposed preconditioning 
technique is tested on a mildly heterogeneous problem 

i7u = - (A + k{x, y)) u = g{x, y). 

The rectangular domain (Om, 600to) x (Om, 1000m) is split in 3 regions where the speed 
of sound c takes 3 different values (1500 m/s, 2000 m/s and 3000 m/s). This brings 
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along heterogeneity in the wave number that is defined as k{x,y) — {2tt f / c{x , y)) , 
with / e {lOHz, 50Hz) the frequency of the point source 



l{x,y) = S{x,y) 



1, if a; = 300, y = 0, 
0, elsewhere, 



The wedge problem was tested for Sommerfeld boundary conditions on all edges of 
the domain with the complex shifted Laplacian preconditioner [5]. For this paper the 
ECS boundary layers are used to absorb outgoing waves, with a complex stretched 
grid preconditioner. 





/ 


10 


20 


30 


40 


50 




interior grid 


32 X 64 


64 X 128 


64 X 128 


128 X 256 


128 X 256 


MG{V{1,0)) 


FGMRES 


38(1.03) 


68(2.57) 


108(5.43) 


128(9.31) 


153(12.74) 




Bi-CGSTAB 


29(0.98) 


64.5(2.17) 


124(4.19) 


134.5(4.64) 


259(8.60) 


MG{V{1,1)) 


FGMRES 


26(0.94) 


44(2.00) 


68(3.56) 


79(5.15) 


96(6.95) 




Bi-CGSTAB 


18(1.04) 


35(2.15) 


53.5(3.29) 


79.5(4.89) 


102(6.30) 



Table 5.2: Number of FGMRES and Bi-CGSTAB iterations (and CPU time per grid 
point xlO"^) for different preconditioning strategies for Case 2. 



In Table [O] we see that MG(V(1,1)) beats MG(V(1,0)) when it comes to outer 
Krylov iterations. The CPU time of the MG(V(1,0)) test is only slightly faster for 
k = 10 and fc = 40 for Bi-CGSTAB. Other tests involved more smoothing steps and 
the use of ILU for approximate inversion as in Case 1, but were again not competitive 
and are therefore left out of the table. 

Case 3: Ionization problem. The model for the last experiment is the 2D Helmholtz 
problem 

Hu = ~{A + k{x, y)) u = g{x, y), 
on the square domain (0, 50)^ with a strongly varying wave number 



1 1 

— H 2 

where 0<A:o<5,0</Lt<10 and a source function 

-f ev 

The south and the west edges of the domain have homogeneous Dirichlet boundary 
conditions while ECS layers absorb outgoing waves on the north and the east edges. 
The problem appears in the simulation of Schrodinger's equation for single and multiple 
ionization of atoms and molecules [20] . The extension to higher dimensions than 2D 
results in a massive amount of storage and computational complexity and is the 
main motivation to develop a robust iterative method for space-dependent Helmholtz 
equations. It was tested for the complex shifted Laplacian and complex stretched grid 
preconditioners in 115]. This model problem is challenging from an iterative point of 
view due the highly space-dependent wave number and the parameter /.t. For values 
^ > 2.73 evanescent waves form near the Dirichlet edges associated to resonances. 
These resonances also appear in the spectrum of the discretized operator as isolated 
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0.5 


1.5 


2.5 


3.5 


4.5 




interior grid 


64^ 


128^ 


256^ 


512^ 


512^ 


MG(V(1,0)) 


FGMRES 


39(0.07) 


127(0.72) 


181(1.54) 


238(3.19) 


281(4.22) 




Bi-CGSTAB 


30.5(0.06) 


303(0.98) 


178(0.62) 


470(1.50) 


565(1.82) 


MG(V(1,1)) 


FGMRES 


24(0.05) 


73(0.42) 


116(0.42) 


137(1.35) 


194(2.40) 




Bi-CGSTAB 


16(0.10) 


62(0.76) 


103(1.08) 


146(1.01) 


227.5(1.57) 



Table 5.3: Number of FGMRES and Bi-CGSTAB iterations (and CPU time per grid 
point xlO^) for different preconditioning strategies for Case 3 with fi = 1. 





ko 


0.5 


1.5 


2.5 


3.5 


4.5 




interior grid 


64^ 


128^ 


256^ 


512^ 


512^ 


MG(V(1,0)) 


FGMRES 


139(0.45) 


159(1.67) 


206(2.71) 


250(3.74) 


285(4.50) 




Bi-CGSTAB 


457.5(0.76) 


294(0.93) 


313(1.05) 


310(1.00) 


562(1.81) 


MG(V(1,1)) 


FGMRES 


75(0.46) 


106(1.52) 


134(1.67) 


150(1.64) 


204(2.73) 




Bi-CGSTAB 


872(3.70) 


138(1.37) 


136.5(1.15) 


186(1.12) 


253(1.43) 



Table 5.4: Number of FGMRES and Bi-CGSTAB iterations (and CPU time per grid 
point xlO^) for different preconditioning strategies for Case 3 with /i = 7. 



eigenvalues on the real axis, left of the smoothest eigenvalues, but are less sensitive 
to grid changes such as complex stretching or grid coarsening. This can pollute the 
multigrid functionality for the preconditioning matrix. 



In Table[5Jwe see that for ko = 1.5 the Bi-CGSTAB method with MG(V(1,0)) 
inversion of the preconditioner takes 303 iterations to converge which is high in 
comparison with FGMRES which only needs 127 iteration. This outlier in the 



convergence is even more apparent when ij, — 7 in Table 5.4 For ko — 0.5 the Bi- 
CGSTAB method with MG(V(1,0)) and MG(V(1,1)) takes 457.5 and 872 iterations 
respectively, in comparison to 139 and 75 for FGMRES that seems to cope better with 
the issues of resonances caused by the coarse grid correction. 

The reason why resonances appear in this problem with a space dependent wave 
number that have evanescent waves, is that at the coarse level the discretization is too 
rough to resolve the details of the variation of the wave number. This leads to slightly 
different eigenvalues on subsequent levels and this gives rise to the resonances in the 
convergences. 

The FMGRES does not suffer from these resonances because the problematic 
eigenmodes are resolved by the Ritz eigenvalue and the Krylov converges superlinearly 
as if this mode is not present. A discussion on superlinear convergence in the iterative 
solution of linear systems can be found in [TO] , 

It is important to note that similar results are obtained when Galerkin is used to 
construct the coarse operators. 

6. Conclusions and outlook. In this paper we have analyzed the iterative 
solution of an Helmholtz equation discretized with finite differences and an absorbing 
boundary condition based on complex scaling. The iterative solver is based on a Krylov 
subspace and it is preconditioned with multigrid that uses a polynomial smoother. 
The multigrid is applied to a Helmholtz problem discretized on a complex valued grid. 

For each level of the multigrid hierarchy, the spectrum of the discrete precondi- 
tioning operator is bounded by a triangle that lies entirely in the lower half of the 
complex plane. Based on the properties of the triangle, we show that it is possible 
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to choose the parameters of the preconditioner such that a third order polynomial 
smoother is stable for all levels and all wave numbers. This smoother can be viewed 
as a sequence of three damped Jacobi steps with three different w's. 

The hand tuned third order polynomial forms an bound for the convergence of 
the GMRES(3) as a smoother. GMRES is more effective, not only because GMRES 
optimizes the coefficients automatically, but also selects the coefficients differently in 
each smoothing step. This will give a better performance than any fixed polynomial 
that is optimized manually. 

Note that GMRES (3) does not necessarily have a smoothing behavior in the strict 
sense. However, the behavior of this GMRES smoother is bounded by the hand tuned 
polynomial that has the smoothing property. We have also found in the various tests 
on model problems that the method gives satisfactory results. 

The numerical results are obtained for the a preconditioner based on complex 
shifted grids, but we expect that similar results will be obtained for a complex shifted 
Laplacian. In jI5j we have shown that the complex shifted grid preconditioner is 
equivalent to the complex shifted Laplacian preconditioner and that it gives the same 
Krylov space convergence. In [5], the complex shifted Laplacian is also inverted with 
ILU. However, we have found that the inversion of the preconditioner based on the 
multigrid combined with the polynomial smoother, performs better, both in compute 
time and required memory than inversion by ILU. 

Unfortunately, the number of Krylov iterations still grows linearly with the wave 
number k and for large wave numbers the number of iterations increases significantly. 



This is illustrated in figure 5.1 We intend to give an analysis of the dependence in an 
upcoming publication. 

We have found that the introduction of complex shifts is not sufficient to eliminate 
all the difficulties with the coarse grid correction. For problems with space dependent 
wave numbers that allow evanescent waves the coarse grid correction can still be 
problematic since the space dependence is felt differently on the various levels of the 
multigrid hierarchy. 

When the multigrid is combined with FGMRES we observe that the convergence 
is not disturbed by diverging coarse grid corrections. This suggests that defiation can 
be effective. This is a possible subject for future research. 

Another subject of future research is the efficient implementation for 3D problems. 
On modern hardware stencil computations are typically communication bound since 
there are only a few floating point operations for each read from the slow memory. 
For the polynomial smoother, however, there are optimizations possible that increase 
the number of floating points operations per load. The application of the polynomial 
smoother can be implemented by loading a cube of the domain into fast memory 
and then applying the three matrix- vector products [6 . In this way we reduce the 
communication with the slow memory. 

Acknowledgements. We acknowledge support from FWO-Flanders through 
grant G. 0174. 08. WV is also supported by FWO-Flanders Krediet aan navorser 
1.5.145.10 N. 



REFERENCES 



[1] A. Bayliss, C.I. Goldstein, and E. Turkel. An iterative method for the Helmholtz equation. 

Journal of Computational Physics, 49(3):443-457, 1983. 
[2] A. Bcn-Mcnahcm and S.J. Singh. Seismic waves and sources. Dover Publications, 1998. 



17 



[3] H. bin Zubair, B. Reps, and W. Vanroose. A preconditioned iterative solver for the scattering 

solutions of the schrodinger equation, accepted for publication in Communications in 

Computational Physics, 2010. 
[4] M. Bollhofcr, M.J. Grotc, and O. Schcnk. Algebraic multilevel preconditioner for the Helmlioltz 

equation in heterogeneous media. SIAM Journal on Scientific Computing, 31(5):3781-3805, 

2009. 

[5] W.L. Briggs, V.E. Hcnson, and S.F. McCormick. A multigrid tutorial. SIAM, 2000. 

[6] K. Datta, S. Kamil, S. Williams, L. Oliker, J. Shall, and K. Yelick. Optimization and performance 

modeling of stencil computations on modern microprocessors. SIAM review, 51(1):129-159, 

2009. 

[7] H.C. Elman, O.G. Ernst, and D.P. O'Leary. A multigrid method enhanced by Krylov sub- 
space iteration for discrete Helmholtz equations. SIAM Journal of Scientific Computing, 
23(4):1291-1315, 2002. 

[8] Y.A. Erlangga, C.W. Oosterlee, and C. Vuik. A novel multigrid based preconditioner for 
heterogeneous Helmholtz problems. SIAM Journal on Scientific Computing, 27(4):1471- 
1492, 2006. 

[9] Y.A. Erlangga, C. Vuik, and C.W. Oosterlee. On a class of preconditioners for solving the 
Helmholtz equation. Applied Numerical Mathematics, 50(3-4):409-425, 2004. 

[10] M.J. Gander and F. Nataf. AILU for Helmholtz problems: a new preconditioner based on the 
analytic parabolic factorization. Journal of Computational Acoustics, 9(4):1499-1506, 2001. 

[11] E. Haber and S. MacLachlan. A fast method for the solution of the helmholtz equation, submitted, 
2010. 

[12] F. Natterer and F. Wubbeling. A propagation-backpropagation method for ultrasound tomogra- 
phy. Inverse Problems, 11:1225, 1995. 

[13] X. Pinel. A peiiAirbcd two-level preeonditioner for the solution of three-dimensional heterogeneous 
Helmlioltz problems witfi applications to geophysics. PhD Thesis, Universitc de Toulouse, 
2010. 

[14] R.E. Plessix and W.A. Mulder. Separation-of-variables as a preconditioner for an iterative 

Helmholtz solver. Applied Numerical Mathematics, 44(3) :385 400. 2003. 
[15] B. Reps, W. Vanroose, and H. bin Zubair. On the indefinite Helmholtz equation: Complex 

stretched absorbing boundary layers, iterative analysis, and preconditioning. Journal of 

Computational Physics, 229(22) :8384-8405, 2010. 
[16] Y. Saad. A flexible inner-outer preconditioned GMRES algorithm. SIAM Journal on Scientific 

Computing, 14:461-461, 1993. 
[17] Y. Shapira. Matrix-Based Multigrid: Theory and Applications. Kluwer Academic Publishers, 

2003. 

[18] H.P. Urbach and D.A. Bernard. Modeling latent-image formation in photolithography, using the 
Helmholtz equation. Journal of the Optical Society of America A, 6(9):1343-1356, 1989. 

[19] H.A. Van der Vorst. Iterative Krylov methods for large linear systems. Cambridge University 
Press, 2003. 

[20] W. Vanroose, F. Martin, T.N. Rescigno, and C.W. McCurdy. Complete photo-induced breakup 
of the h2 molecule as a probe of molecular electron correlation. Science, 310(5755):1787, 
2005. 



18 



0.5 




1 ■ 10^ 



0.25 



-0.25 



-0.5 



0.02 



0.018 



0.016 



0.014 



0.012 





0.9994 



0.9996 



0.9998 
5ft/(A) 



Fig. 3.1: This figure illustrates Example [T] a) The spectrum of the Helmholtz problem 
on a complex stretched grid for k — 100 on the domain (0,1)'^. It is discretized 
with 192 by 192 finite difference points on the interior and ECS with 64 points as an 
absorbing boundary layer. The grid is rotated over 8 degrees in the interior and 30 
degrees on the exterior. The spectrum is bounded by the line segments that connect 
— fc^, and Ajj*^-* and A^'^'' that form a triangle, b) Map of the spectrum under the third 

order polynomial f{t) which is constructed such that f{~k'^) = e"^, ,f{X^p^) = and 

/(A^^'=^) = 0. For this map if — 0.9167 degrees, c) The smoothest modes are mapped 
near the unit circle. Depending on the angle Lp, these are mapped outside or inside 
the unit circle. 
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Fig. 3.2: Illustration of the conditions in (3.11) and (3.161 in a multigrid hierarchy 
and as a function of the wave number k and the rotation angle 6fj of the mesh width 
of the interior region. The details of the multigrid hierarchy is given in Example |2] 
The dashed lines are drawn where |5ft(62/^|)| = |^i/^^| for each of the levels. For each 
level and wave numbers on the right of the dashed line, any choice of the rotation 
angle Op and gives a stable smoother for that particular level. On the left of the 
dashed line we need to solve for tp. However, only the region above the solid line gives 
a stable smoother. The solid lines are based on the condition in (3.161. In order to 
build a stable smoother on each level of the hierarchy for all possible k we can choose 
the angle of rotation of the inner region 9p around 10 degrees, which is above all solid 
lines. 
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Fig. 3.3: This figure illustrates that it is not sufRcient to choose the angle Lp such 
that the top of the bounding triangle of the spectrum is tangent to the unit circle. 
Although the smoothest modes are still mapped inside the circle and the oscillatory 
modes are mapped near zero, an intermediate part of the spectrum goes outside the 
circle under the transformation, leading to unstable modes. The solution is to adjust 
/i^, the mesh width in the interior region, until the complete spectrum is mapped 
inside the circle. In Lemma 13.31 a condition is formulated that realizes this. The solid 
lines form the image of the bounding triangle. The 2D example is discretized with 32 
grid points in the interior and 16 grid points in the ECS layer in one direction. The 
interior is scaled by an angle of 5 degrees, while the exterior has 30 degrees. The wave 
number is k=40. 
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Fig. 3.4: The map of the spectrum under the third order polynomial smoother for 
each level of the multigrid hierarchy. The wave number is fc = 40. The interior angle 
of the complex stretched grid is 0^ = 7r/18 , which leads to a stable smoother for all 
levels and all k. The ECS rotation is 9-^ ~ tt/G. The grid is a Kronecker product of a 
grid with 7, 15, 31, 63, 127 and 255 points in one dimension. One quarter of the grid 
points lie on the ECS contour. 
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Fig. 4.1: The measured two grid convergence with one pre and one post smoothing 
step for the Helmholtz equation on a (0, 1)^ domain as a function of the wave number k. 
The domain has a ECS layer on two of the boundaries. We compare the convergence 
of the polynomial smoother with the convergence of GMRES(l), GMRES(2) and 
GMRES(3) as smoother. From the top left to the bottom right we show results for 
31^, 63^, 127^, 255^ and 511^ grid points. The rotation angle of the inner part is 
0^ = Yg = 0.1745, which is above the critical angle of the polynomial smoother. We 
note that for each of these two grid tests the convergence suffers when k is such that 

(k) 

X'^"-' of the coarser level, crosses the imaginary axis. Since the polynomial smoother is 
a particular polynomial, it gives a worst case bound for the GMRES(3) smoother. 
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Fig. 4.2: The measured convergence rate of different V-cyles with GMRES(3)- 
smoothing apphed to the preconditioning matrix M^^'^ , as a function of the wave 
number k. The matrix is buih by discretizing the homogeneous Helmhohz problem 
with ECS boundary layers on a complex stretched grid with Op — 0.18 as the rotation 
angle for the interior domain. The dashed lines show the residual reduction after one 
V-cycle, while the solid lines are average rates over more subsequent cycles. Doubling 
the amount of smoothing from V(1,0) (+) to V(l,l) (*) improves the rate from under 
40% to under 25%. The performance of the first V-cycle alone even decreases from un- 
der 35% to under 8%. Eventually only one such V-cycle will be used to approximately 
invert the preconditioner Mj^^'~^. The experiment with V(2,l) (x) indicates that an 
extra smoothing step does not seem to improve the rates in the same order, maximum 
15% average and maximum 5% for the first cycle. 
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Fig. 5.1: Number of iterations (left) and CPU time scaled by the total number of 
grid points (right) for FGMRES and Bi-CGSTAB as a function of the wave number k. 
The Helmholtz matrix H was preconditioned with one V(l,l) cycle with GMRES(3)- 
smoothing applied to the M^^'^ with angle Op = 0.18. The number of grid points 
grows with k according to the rule kh < 0.625. 
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